#!/usr/bin/env python
# coding: utf-8

# # Indus_5M

# In[1]:


import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import math
#from scipy.stats import skew


# In[2]:


df = pd.read_csv("D:/research/recordbreak/model_streamflow/indus45.csv")


# In[3]:


print(df)


# In[4]:


df.set_index('date', inplace=True)
df.index = pd.to_datetime(df.index)


# In[5]:


df['year'] = df.index.year
start_date = '1962-01-01'
end_date = '2099-12-31'
df = df[start_date:end_date]
grouped = df.groupby('year')


# In[6]:


print(df)


# In[7]:


# Select the 'obs' and 'water year' columns
#df_to_save = df[['obs', 'water year']]

# Save the selected columns to a CSV file with the index
#df_to_save.to_csv('obs_and_water_year.csv', index=True)


# In[8]:


#obs_max.to_csv('obs_max.csv', index=True)


# # Obs

# In[9]:


obs_max = grouped['obs'].max()
print(obs_max)


# In[10]:


# Initialize an empty list to store the binary values
obs_bin = []

# Initialize the current maximum value to a very small number
current_max = float('-inf')

# Initialize a counter
counter = 0

# Iterate through the maximum values
for value in obs_max:
    if value > current_max:
        current_max = value
        counter += 1
    obs_bin.append(counter)

# Convert the list to a Pandas Series
obs_bin = pd.Series(obs_bin, index=obs_max.index)

# Print the binary series
print(obs_bin)


# In[11]:


# Assuming 'n' is defined
n= [i+1 for i in range(len(obs_bin))]  # Generate a list of n values

# fr=flood record breaking events decay rate
fr_obs = []

# Iterate through the obs_bin values and corresponding n values
for obs_bin, n in zip(obs_bin, n):
    fr_obs.append(obs_bin / n)

# Convert the list to a Pandas Series
fr_obs = pd.Series(fr_obs, index=obs_max.index)

# Print the resulting divided series
print(fr_obs)


# In[12]:


# Generate n values from 1 to 61 ....139
n_values = list(range(1,139))

# Calculate the theoretical flood record-breaking events decay rate
fr_theo = 1 / pd.Series(n_values, index=obs_max.index)

# Print the resulting theoretical series
print(fr_theo)


# In[13]:


fr_obs.index = fr_obs.index.astype(str)

plt.figure(figsize=(10, 6))

# Plotting the flood record-breaking events decay rate
plt.plot(fr_obs.index, fr_obs.values, label='Flood Record-breaking Events Decay Rate', color='green')
plt.plot( fr_theo.values, label='theoratical line', color='blue',linestyle='--')
# Adding labels and title
plt.xlabel('Water Year')
plt.ylabel('Decay Rate')
plt.title('Flood Record-breaking Events Decay Rate Over Time')
plt.xticks(fr_obs.index[::4])
#plt.ylim(-0.5, 1.5)
# Adding a legend
plt.legend()

# Display the plot
plt.show()


# # Sim1

# In[14]:


sim1_max=grouped['sim1'].max()
print(sim1_max)


# In[15]:


#dnn_max_df = dnn_max.to_frame(name='Max_Values')

# Save the DataFrame to a CSV file
#dnn_max_df.to_csv('dnn_max_values.csv', index=True)


# In[16]:


# Initialize an empty list to store the binary values
sim1_bin = []

# Initialize the current maximum value to a very small number
current_max = float('-inf')

# Initialize a counter
counter = 0

# Iterate through the maximum values
for value in sim1_max:
    if value > current_max:
        current_max = value
        counter += 1
    sim1_bin.append(counter)

# Convert the list to a Pandas Series
sim1_bin = pd.Series(sim1_bin, index=sim1_max.index)

# Print the resulting binary series
print(sim1_bin)


# In[17]:


#sim1_bin.to_csv('sim1_bin.csv', index=True)


# In[18]:


# Assuming 'n' is defined
n= [i+1 for i in range(len(sim1_bin))]  # Generate a list of n values

# fr=flood record breaking events decay rate
fr_sim1 = []

# Iterate through the obs_bin values and corresponding n values
for sim1_bin, n in zip(sim1_bin, n):
    fr_sim1.append(sim1_bin / n)

# Convert the list to a Pandas Series
fr_sim1 = pd.Series(fr_sim1, index=sim1_max.index)

# Print the resulting divided series
print(fr_sim1)


# # sim2

# In[19]:


sim2_max=grouped['sim2'].max()
print(sim2_max)


# In[20]:


# Initialize an empty list to store the binary values
sim2_bin = []

# Initialize a counter
counter = 0

# Initialize the current maximum value to a very small number
current_max = float('-inf')


# Iterate through the maximum values
for value in sim2_max:
    if value > current_max:
        current_max = value
        counter += 1
    sim2_bin.append(counter)

# Convert the list to a Pandas Series
sim2_bin = pd.Series(sim2_bin, index=sim2_max.index)

# Print the resulting binary series
print(sim2_bin)


# In[21]:


# Assuming 'n' is defined
n= [i+1 for i in range(len(sim2_bin))]  # Generate a list of n values

# fr=flood record breaking events decay rate
fr_sim2 = []

# Iterate through the obs_bin values and corresponding n values
for sim2_bin, n in zip(sim2_bin, n):
    fr_sim2.append(sim2_bin / n)

# Convert the list to a Pandas Series
fr_sim2 = pd.Series(fr_sim2, index=sim2_max.index)

# Print the resulting divided series
print(fr_sim2)


# # Sim3

# In[22]:


sim3_max=grouped['sim3'].max()
print(sim3_max)


# In[23]:


# Initialize an empty list to store the binary values
sim3_bin = []

# Initialize a counter
counter = 0

# Initialize the current maximum value to a very small number
current_max = float('-inf')


# Iterate through the maximum values
for value in sim3_max:
    if value > current_max:
        current_max = value
        counter += 1
    sim3_bin.append(counter)

# Convert the list to a Pandas Series
sim3_bin = pd.Series(sim3_bin, index=sim3_max.index)

# Print the resulting binary series
print(sim3_bin)


# In[24]:


# Assuming 'n' is defined
n= [i+1 for i in range(len(sim3_bin))]  # Generate a list of n values

# fr=flood record breaking events decay rate
fr_sim3 = []

# Iterate through the obs_bin values and corresponding n values
for sim3_bin, n in zip(sim3_bin, n):
    fr_sim3.append(sim3_bin / n)

# Convert the list to a Pandas Series
fr_sim3 = pd.Series(fr_sim3, index=sim3_max.index)

# Print the resulting divided series
print(fr_sim3)


# # sim4

# In[25]:


sim4_max=grouped['sim4'].max()
print(sim4_max)


# In[26]:


# Initialize an empty list to store the binary values
sim4_bin = []

# Initialize a counter
counter = 0

# Initialize the current maximum value to a very small number
current_max = float('-inf')


# Iterate through the maximum values
for value in sim4_max:
    if value > current_max:
        current_max = value
        counter += 1
    sim4_bin.append(counter)

# Convert the list to a Pandas Series
sim4_bin = pd.Series(sim4_bin, index=sim4_max.index)

# Print the resulting binary series
print(sim4_bin)


# In[27]:


# Assuming 'n' is defined
n= [i+1 for i in range(len(sim4_bin))]  # Generate a list of n values

# fr=flood record breaking events decay rate
fr_sim4 = []

# Iterate through the obs_bin values and corresponding n values
for sim4_bin, n in zip(sim4_bin, n):
    fr_sim4.append(sim4_bin / n)

# Convert the list to a Pandas Series
fr_sim4 = pd.Series(fr_sim4, index=sim4_max.index)

# Print the resulting divided series
print(fr_sim4)


# # Sim5

# In[28]:


sim5_max=grouped['sim5'].max()
print(sim5_max)


# In[29]:


# Initialize an empty list to store the binary values
sim5_bin = []

# Initialize a counter
counter = 0

# Initialize the current maximum value to a very small number
current_max = float('-inf')


# Iterate through the maximum values
for value in sim5_max:
    if value > current_max:
        current_max = value
        counter += 1
    sim5_bin.append(counter)

# Convert the list to a Pandas Series
sim5_bin = pd.Series(sim5_bin, index=sim5_max.index)

# Print the resulting binary series
print(sim5_bin)


# In[30]:


# Assuming 'n' is defined
n= [i+1 for i in range(len(sim5_bin))]  # Generate a list of n values

# fr=flood record breaking events decay rate
fr_sim5 = []

# Iterate through the obs_bin values and corresponding n values
for sim5_bin, n in zip(sim5_bin, n):
    fr_sim5.append(sim5_bin / n)

# Convert the list to a Pandas Series
fr_sim5 = pd.Series(fr_sim5, index=sim5_max.index)

# Print the resulting divided series
print(fr_sim5)


# In[ ]:





# # Plot

# In[31]:


plt.figure(figsize=(20, 10))

# Plotting the flood record-breaking events decay rate
#plt.plot(fr_obs.index, fr_obs.values, label='Observed', color='black',alpha=0.8,linewidth=2.5)
plt.plot(fr_obs['1962':'2022'].index, fr_obs['1962':'2022'].values, label='Observed', color='black', alpha=0.8, linewidth=5)
plt.plot(fr_obs.index,fr_sim1.values,label='ESM2',linestyle='dotted',color='red',linewidth=4)
plt.plot(fr_obs.index,fr_sim2.values,  label='MK3.6',linestyle='dashdot',color='magenta',linewidth=4)
plt.plot(fr_obs.index,fr_sim3.values, label='ESM2M',linestyle='dashed',color='brown',linewidth=4)
plt.plot(fr_obs.index,fr_sim4.values,label='CM5A',linestyle='dotted',color='green',linewidth=4)
plt.plot(fr_obs.index,fr_sim5.values, label='ESM',linestyle='dashdot',color='orange',linewidth=4)
plt.plot( fr_obs.index,fr_theo.values, label='Theoratical(1/n)', color='blue',linewidth=4)
# Adding labels and title
#plt.xlabel('Water Year (Apr-Mar)',fontsize=16,fontweight='bold')
#plt.ylabel('Probability (log)',fontsize=18)
#plt.title('Indus wetting Annual-Max Record-breaking',fontsize=20,fontweight='bold')
plt.xticks(fr_obs.index[::17])

plt.yscale('log')
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02,0.01]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)


plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.gca().tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=16, width=4)
#plt.ylim(-0.5, 1.5)
# Adding a legend
plt.legend(loc='upper right',prop={'size': 22})
#plt.xlim(1960, 2099)
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)

# Set x-axis limits within the data range
#plt.xlim(min(fr_obs.index), max(fr_obs.index))
# Display the plot
plt.show()


# In[32]:


plt.figure(figsize=(20, 10))

# Plotting the flood record-breaking events decay rate
#plt.plot(fr_obs.index, fr_obs.values, label='Observed', color='black',alpha=0.8,linewidth=2.5)
plt.plot(fr_obs['1962':'2022'].index, fr_obs['1962':'2022'].values, label='Observed', color='black', alpha=0.8, linewidth=5)
plt.plot(fr_obs.index,fr_sim1.values,label='CanESM2',linestyle='dotted',color='red',linewidth=4)
plt.plot(fr_obs.index,fr_sim2.values,  label='CSIRO-MK3.6',linestyle='dashdot',color='magenta',linewidth=4)
plt.plot(fr_obs.index,fr_sim3.values, label='GFDL-ESM2M',linestyle='dashed',color='brown',linewidth=4)
plt.plot(fr_obs.index,fr_sim4.values,label='IPSL-CM5A-LR',linestyle='dotted',color='green',linewidth=4)
plt.plot(fr_obs.index,fr_sim5.values, label='MPI-ESM-MR',linestyle='dashdot',color='orange',linewidth=4)
plt.plot( fr_obs.index,fr_theo.values, label='Theoratical (1/n)', color='blue',linewidth=3)
# Adding labels and title
#plt.xlabel('Water Year (Apr-Mar)',fontsize=16,fontweight='bold')
#plt.ylabel('Probability (log)',fontsize=18)
#plt.title('Indus wetting Annual-Max Record-breaking',fontsize=20,fontweight='bold')
plt.xticks(fr_obs.index[::17])

plt.yscale('log')
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)


plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.gca().tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=16, width=4)
#plt.ylim(-0.5, 1.5)
# Adding a legend
#legend = plt.legend(loc='upper right', prop={'size': 22}, bbox_to_anchor=(0.99, 0.98), borderaxespad=0.)
#legend.get_frame().set_linewidth(4)
#legend.get_frame().set_edgecolor('black')

ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)

# Set x-axis limits within the data range
plt.xlim(min(fr_obs.index), max(fr_obs.index))
# Display the plot
#plt.xlim(1961, 2022)
plt.show()


# # Min Record Breaking 

# In[33]:


obs_min = grouped['obs'].min()
print(obs_min)


# In[34]:


#obs_min.to_csv('obs_min.csv', index=True)


# In[35]:


# Initialize an empty list to store the binary values
obs_min_bin = []

# Initialize the current minimum value to a very large number
current_min = float('-inf')

# Initialize a counter
counter = 0

# Iterate through the minimum values
for value in obs_min:
    if value > current_min:
        current_min = value
        counter += 1
    obs_min_bin.append(counter)

# Convert the list to a Pandas Series
obs_min_bin = pd.Series(obs_min_bin, index=obs_min.index)

# Print the resulting binary series
print(obs_min_bin)


# In[36]:


#obs_min_bin.to_csv('obs_min_bin.csv', index=True)


# In[37]:


# Assuming 'n' is defined
n= [i+1 for i in range(len(obs_min_bin))]  # Generate a list of n values

# fr=flood record breaking events decay rate
dr_obs = []

# Iterate through the obs_bin values and corresponding n values
for obs_min_bin, n in zip(obs_min_bin, n):
    dr_obs.append(obs_min_bin / n)

# Convert the list to a Pandas Series
dr_obs = pd.Series(dr_obs, index=obs_min.index)

# Print the resulting divided series
print(dr_obs)


# In[39]:


# Generate n values from 1 to 61
n_values = list(range(1, 139))

# Calculate the theoretical flood record-breaking events decay rate
dr_theo = 1 / pd.Series(n_values, index=obs_min.index)

# Print the resulting theoretical series
print(dr_theo)


# In[40]:


dr_obs.index = dr_obs.index.astype(str)

plt.figure(figsize=(10, 6))

# Plotting the flood record-breaking events decay rate
plt.plot(dr_obs.index, dr_obs.values, label='Flood Record-breaking Events Decay Rate', color='green')
plt.plot( dr_theo.values, label='theoratical line', color='blue',linestyle='--')
# Adding labels and title
plt.xlabel('Water Year')
plt.ylabel('Decay Rate')
plt.title('Flood Record-breaking Events Decay Rate Over Time')
plt.xticks(fr_obs.index[::4])
#plt.ylim(-0.5, 1.5)
# Adding a legend
plt.legend()

# Display the plot
plt.show()


# # Sim1

# In[41]:


sim1_min=grouped['sim1'].min()
print(sim1_min)


# In[42]:


#dnn_min_df = dnn_min.to_frame(name='min_Values')

# Save the DataFrame to a CSV file
#dnn_min_df.to_csv('dnn_min_values.csv', index=True)


# In[43]:


# Initialize an empty list to store the binary values
sim1_min_bin = []


# Initialize the current minimum value to a very large number
current_min = float('-inf')

# Initialize a counter
counter = 0

# Iterate through the minimum values
for value in sim1_min:
    if value > current_min:
        current_min = value
        counter += 1
    sim1_min_bin.append(counter)

# Convert the list to a Pandas Series
sim1_min_bin = pd.Series(sim1_min_bin, index=sim1_min.index)

# Print the resulting binary series
print(sim1_min_bin)


# In[44]:


# Assuming 'n' is defined
n= [i+1 for i in range(len(sim1_min_bin))]  # Generate a list of n values

# fr=flood record breaking events decay rate
dr_sim1 = []

# Iterate through the obs_bin values and corresponding n values
for sim1_min_bin, n in zip(sim1_min_bin, n):
    dr_sim1.append(sim1_min_bin / n)

# Convert the list to a Pandas Series
dr_sim1 = pd.Series(dr_sim1, index=sim1_min.index)

# Print the resulting divided series
print(dr_sim1)


# # Sim2

# In[45]:


sim2_min=grouped['sim2'].min()
print(sim2_min)


# In[46]:


# Initialize an empty list to store the binary values
sim2_min_bin = []

current_min = float('-inf')

# Initialize a counter
counter = 0

# Iterate through the minimum values
for value in sim2_min:
    if value > current_min:
        current_min = value
        counter += 1
    sim2_min_bin.append(counter)

# Convert the list to a Pandas Series
sim2_min_bin = pd.Series(sim2_min_bin, index=sim2_min.index)

# Print the resulting binary series
print(sim2_min_bin)


# In[47]:


# Assuming 'n' is defined
n= [i+1 for i in range(len(sim2_min_bin))]  # Generate a list of n values

# fr=flood record breaking events decay rate
dr_sim2 = []

# Iterate through the obs_bin values and corresponding n values
for sim2_min_bin, n in zip(sim2_min_bin, n):
    dr_sim2.append(sim2_min_bin / n)

# Convert the list to a Pandas Series
dr_sim2 = pd.Series(dr_sim2, index=sim2_min.index)

# Print the resulting divided series
print(dr_sim2)


# # Sim3

# In[48]:


sim3_min=grouped['sim3'].min()
print(sim3_min)


# In[49]:


# Initialize an empty list to store the binary values
sim3_min_bin = []

current_min = float('-inf')

# Initialize a counter
counter = 0

# Iterate through the minimum values
for value in sim3_min:
    if value > current_min:
        current_min = value
        counter += 1
    sim3_min_bin.append(counter)

# Convert the list to a Pandas Series
sim3_min_bin = pd.Series(sim3_min_bin, index=sim3_min.index)

# Print the resulting binary series
print(sim3_min_bin)


# In[50]:


# Assuming 'n' is defined
n= [i+1 for i in range(len(sim3_min_bin))]  # Generate a list of n values

# fr=flood record breaking events decay rate
dr_sim3 = []

# Iterate through the obs_bin values and corresponding n values
for sim3_min_bin, n in zip(sim3_min_bin, n):
    dr_sim3.append(sim3_min_bin / n)

# Convert the list to a Pandas Series
dr_sim3 = pd.Series(dr_sim3, index=sim3_min.index)

# Print the resulting divided series
print(dr_sim3)


# # Sim4

# In[51]:


sim4_min=grouped['sim4'].min()
print(sim4_min)


# In[52]:


# Initialize an empty list to store the binary values
sim4_min_bin = []

current_min = float('-inf')

# Initialize a counter
counter = 0

# Iterate through the minimum values
for value in sim4_min:
    if value > current_min:
        current_min = value
        counter += 1
    sim4_min_bin.append(counter)

# Convert the list to a Pandas Series
sim4_min_bin = pd.Series(sim4_min_bin, index=sim4_min.index)

# Print the resulting binary series
print(sim4_min_bin)


# In[53]:


# Assuming 'n' is defined
n= [i+1 for i in range(len(sim4_min_bin))]  # Generate a list of n values

# fr=flood record breaking events decay rate
dr_sim4 = []

# Iterate through the obs_bin values and corresponding n values
for sim4_min_bin, n in zip(sim4_min_bin, n):
    dr_sim4.append(sim4_min_bin / n)

# Convert the list to a Pandas Series
dr_sim4 = pd.Series(dr_sim4, index=sim4_min.index)

# Print the resulting divided series
print(dr_sim4)


# # sim5

# In[54]:


sim5_min=grouped['sim5'].min()
print(sim5_min)


# In[55]:


# Initialize an empty list to store the binary values
sim5_min_bin = []

# Initialize a counter
current_min = float('-inf')

# Initialize a counter
counter = 0

# Iterate through the minimum values
for value in sim5_min:
    if value > current_min:
        current_min = value
        counter += 1
    sim5_min_bin.append(counter)

# Convert the list to a Pandas Series
sim5_min_bin = pd.Series(sim5_min_bin, index=sim5_min.index)

# Print the resulting binary series
print(sim5_min_bin)


# In[56]:


# Assuming 'n' is defined
n= [i+1 for i in range(len(sim5_min_bin))]  # Generate a list of n values

# fr=flood record breaking events decay rate
dr_sim5 = []

# Iterate through the obs_bin values and corresponding n values
for sim5_min_bin, n in zip(sim5_min_bin, n):
    dr_sim5.append(sim5_min_bin / n)

# Convert the list to a Pandas Series
dr_sim5 = pd.Series(dr_sim5, index=sim5_min.index)

# Print the resulting divided series
print(dr_sim5)


# # Plot

# In[57]:


plt.figure(figsize=(20, 10))

# Plotting the flood record-breaking events decay rate
#plt.plot(dr_obs.index, dr_obs.values, label='Observed', color='black',alpha=0.9,linewidth=2.5)
plt.plot(dr_obs['1962':'2022'].index, dr_obs['1962':'2022'].values, label='Observed', color='black', alpha=0.9, linewidth=5)
plt.plot(dr_obs.index,dr_sim1.values,label='ESM2',linestyle='dotted',color='red',linewidth=4)
plt.plot(dr_obs.index,dr_sim2.values,  label='MK3.6',linestyle='dashdot',color='magenta',linewidth=4)
plt.plot(dr_obs.index,dr_sim3.values, label='ESM2M',linestyle='dashed',color='brown',linewidth=4)
plt.plot(dr_obs.index,dr_sim4.values,label='CM5A',linestyle='dotted',color='green',linewidth=4)
plt.plot(dr_obs.index,dr_sim5.values, label='ESM',linestyle='dashdot',color='orange',linewidth=4)
plt.plot( dr_obs.index,dr_theo.values, label='Theo (1/n)', color='blue',linewidth=3)
# Adding labels and title
#plt.xlabel('Water Year (Apr-Mar)',fontsize=16,fontweight='bold')
#plt.ylabel('Record Rate (log)',fontsize=18)
#plt.title('Indus wetting Annual-Min Record-breaking',fontsize=20,fontweight='bold')
plt.xticks(fr_obs.index[::17])

plt.yscale('log')
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)


plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.gca().tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=16, width=4)
#plt.ylim(-0.5, 1.5)
# Adding a legend
#legend = plt.legend(loc='upper right', prop={'size': 26}, bbox_to_anchor=(0.99, 0.98), borderaxespad=0.)
#legend.get_frame().set_linewidth(4)
#legend.get_frame().set_edgecolor('black')
#plt.xlim(1960, 2099)
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)

# Set x-axis limits within the data range
plt.xlim(min(dr_obs.index), max(dr_obs.index))
# Display the plot
plt.show()


# In[58]:


plt.figure(figsize=(20, 10))

# Plotting the flood record-breaking events decay rate
#plt.plot(fr_obs.index, fr_obs.values, label='Observed', color='black',alpha=0.8,linewidth=2.5)
plt.plot(dr_obs['1962':'2022'].index, dr_obs['1962':'2022'].values, label='Observed', color='black', alpha=0.9, linewidth=5)
plt.plot(dr_obs.index,dr_sim1.values,label='CanESM2',linestyle='dotted',color='red',linewidth=4)
plt.plot(dr_obs.index,dr_sim2.values,  label='CSIRO-MK3.6',linestyle='dashdot',color='magenta',linewidth=4)
plt.plot(dr_obs.index,dr_sim3.values, label='GFDL-ESM2M',linestyle='dashed',color='brown',linewidth=4)
plt.plot(dr_obs.index,dr_sim4.values,label='IPSL-CM5A-LR',linestyle='dotted',color='green',linewidth=4)
plt.plot(dr_obs.index,dr_sim5.values, label='MPI-ESM-MR',linestyle='dashdot',color='orange',linewidth=4)
plt.plot( dr_obs.index,dr_theo.values, label='Theoratical (1/n)', color='blue',linewidth=4)
# Adding labels and title
#plt.xlabel('Water Year (Apr-Mar)',fontsize=16,fontweight='bold')
#plt.ylabel('Probability (log)',fontsize=18)
#plt.title('Indus wetting Annual-Min Record-breaking',fontsize=20,fontweight='bold')
plt.xticks(fr_obs.index[::17])

plt.yscale('log')
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02,0.01]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)


plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.gca().tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=16, width=4)
#plt.ylim(-0.5, 1.5)
# Adding a legend
legend = plt.legend(loc='upper right', prop={'size': 26}, bbox_to_anchor=(0.99, 0.98), borderaxespad=0.)
legend.get_frame().set_linewidth(4)
legend.get_frame().set_edgecolor('black')

ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)

# Set x-axis limits within the data range
plt.xlim(min(dr_obs.index), max(dr_obs.index))
# Display the plot
#plt.xlim(1961, 2022)
plt.show()


# # Ratio

# In[59]:


print(fr_obs)


# In[60]:


print(dr_obs)


# In[61]:


# Calculate the ratio between fr_obs and dr_obs
obs_ratio = fr_obs / dr_obs

# Print the resulting ratio series
print(obs_ratio)


# In[62]:


# Calculate the ratio between fr_lstm and dr_lstm
sim2_ratio = fr_sim2 / dr_sim2

# Print the resulting ratio series
print(sim2_ratio)


# In[63]:


# Calculate the ratio between fr_dnn and dr_dnn
sim1_ratio = fr_sim1 / dr_sim1

# Print the resulting ratio series
print(sim1_ratio)


# In[64]:


# Calculate the ratio between fr_gru and dr_gru
sim3_ratio = fr_sim3 / dr_sim3

# Print the resulting ratio series
print(sim3_ratio)


# In[65]:


# Calculate the ratio between fr_gru and dr_gru
sim4_ratio = fr_sim4 / dr_sim4

# Print the resulting ratio series
print(sim4_ratio)


# In[66]:


# Calculate the ratio between fr_gru and dr_gru
sim5_ratio = fr_sim5 / dr_sim5

# Print the resulting ratio series
print(sim5_ratio)


# # Ratio Plot

# In[80]:


plt.figure(figsize=(18, 10))

# Plotting the flood record-breaking events decay rate
#plt.plot(dr_obs.index, obs_ratio.values, label='Observed', color='black',alpha=0.9,linewidth=2.5)
plt.plot(dr_obs['1962':'2022'].index, obs_ratio['1962':'2022'].values, label='Observed', color='black', alpha=0.9, linewidth=5.5)
plt.plot(dr_obs.index,sim1_ratio.values,label='CanESME2',linestyle='dotted',color='red',linewidth=4.5)
plt.plot(dr_obs.index,sim2_ratio.values,  label='CSIRO-MK3.6',linestyle='dashdot',color='magenta',linewidth=4.5)
plt.plot(dr_obs.index,sim3_ratio.values, label='GFDL-ESM2M',linestyle='dashed',color='brown',linewidth=4.5)
plt.plot(dr_obs.index,sim4_ratio.values,label='IPSL-CM5A-LR',linestyle='dotted',color='green',linewidth=4.5)
plt.plot(dr_obs.index,sim5_ratio.values, label='MPI-ESM-MR',linestyle='dashdot',color='orange',linewidth=4.5)
#Theoratical
plt.axhline(y=1, color='blue',label='1/n',  linestyle='solid',linewidth=3)

#plt.plot( fr_theo.values, label='theoratical', color='black',linestyle='--')
# Adding labels and title
#plt.xlabel('Water Year (Apr-Mar)',fontsize=16,fontweight='bold')
#plt.ylabel('Record Ratio (log)',fontsize=16)
#plt.title('Indus wetting Annual Max-Min Record-breaking Ratio', fontsize=20, fontweight='bold')

plt.xticks(fr_obs.index[::15])
plt.yscale('log')

yticks_log = [0.2, 0.5, 1, 2, 5, 10]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)


plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.gca().tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=16, width=4)
#plt.ylim(-0.5, 1.5)
# Adding a legend
legend = plt.legend(loc='upper right', prop={'size': 30}, bbox_to_anchor=(0.99, 0.98), borderaxespad=0.,ncol=2)
legend.get_frame().set_linewidth(4)
legend.get_frame().set_edgecolor('black')
#plt.xlim(1960, 2099)
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)

# Set x-axis limits within the data range
plt.xlim(min(dr_obs.index), max(dr_obs.index))
#plt.xlim(1960, 2100)
# Display the plot
plt.show()


# # WGW_Max

# In[79]:


plt.figure(figsize=(20, 10))

# Plotting the flood record-breaking events decay rate
#plt.plot(fr_obs.index, fr_obs.values, label='Observed', color='black',alpha=0.8,linewidth=2.5)
plt.plot(fr_obs['1962':'2022'].index, fr_obs['1962':'2022'].values, label='Observed', color='black', alpha=0.9, linewidth=5)
plt.plot(fr_obs.index,fr_sim1.values,label='CanESM2',linestyle='dotted',color='red',linewidth=4)
plt.plot(fr_obs.index,fr_sim2.values,  label='CSIRO-MK3.6',linestyle='dashdot',color='magenta',linewidth=4)
plt.plot(fr_obs.index,fr_sim3.values, label='GFDL-ESM2M',linestyle='dashed',color='brown',linewidth=4)
plt.plot(fr_obs.index,fr_sim4.values,label='IPSL-CM5A-LR',linestyle='dotted',color='green',linewidth=4)
plt.plot(fr_obs.index,fr_sim5.values, label='MPI-ESM-MR',linestyle='dashdot',color='orange',linewidth=4)
plt.plot( fr_obs.index,fr_theo.values, label='Theoratical (1/n)', color='blue',linewidth=3, alpha=0.8)
# Adding labels and title
#plt.xlabel('Water Year (Apr-Mar)',fontsize=16,fontweight='bold')
#plt.ylabel('Probability (log)',fontsize=18)
#plt.title('Indus wetting Annual-Max Record-breaking',fontsize=20,fontweight='bold')
plt.xticks(fr_obs.index[::15])

plt.yscale('log')
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02,0.01]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)


plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.gca().tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=16, width=4)
#plt.ylim(-0.5, 1.5)
# Adding a legend
#legend = plt.legend(loc='upper right', prop={'size': 22}, bbox_to_anchor=(0.99, 0.98), borderaxespad=0.)
#legend.get_frame().set_linewidth(4)
#legend.get_frame().set_edgecolor('black')

ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)

# Set x-axis limits within the data range
plt.xlim(min(fr_obs.index), max(fr_obs.index))
# Display the plot
#plt.xlim(1961, 2022)
plt.show()


# # WGW_Min

# In[69]:


plt.figure(figsize=(20, 10))

# Plotting the flood record-breaking events decay rate
#plt.plot(fr_obs.index, fr_obs.values, label='Observed', color='black',alpha=0.8,linewidth=2.5)
plt.plot(dr_obs['1962':'2022'].index, dr_obs['1962':'2022'].values, label='Observed', color='black', alpha=0.9, linewidth=5)
plt.plot(dr_obs.index,dr_sim1.values,label='CanESM2',linestyle='dotted',color='red',linewidth=4)
plt.plot(dr_obs.index,dr_sim2.values,  label='CSIRO-MK3.6',linestyle='dashdot',color='magenta',linewidth=4)
plt.plot(dr_obs.index,dr_sim3.values, label='GFDL-ESM2M',linestyle='dashed',color='brown',linewidth=4)
plt.plot(dr_obs.index,dr_sim4.values,label='IPSL-CM5A-LR',linestyle='dotted',color='green',linewidth=4)
plt.plot(dr_obs.index,dr_sim5.values, label='MPI-ESM-MR',linestyle='dashdot',color='orange',linewidth=4)
plt.plot( dr_obs.index,dr_theo.values, label='Theoratical (1/n)', color='blue',linewidth=3, alpha=0.8)
# Adding labels and title
#plt.xlabel('Water Year (Apr-Mar)',fontsize=16,fontweight='bold')
#plt.ylabel('Probability (log)',fontsize=18)
#plt.title('Indus wetting Annual-Min Record-breaking',fontsize=20,fontweight='bold')
plt.xticks(fr_obs.index[::15])

plt.yscale('log')
yticks_log = [1, 0.5,0.2,0.1,0.05,0.02,0.01]  # Adjust these values based on your desired logarithmic scale
yticklabels = [str(val) for val in yticks_log]  # Convert values to strings for labeling
plt.yticks(yticks_log, yticklabels)


plt.xticks(fontsize=30)
plt.yticks(fontsize=30)
plt.gca().tick_params(axis='x', which='major', length=16, width=4)  # Set the desired length and width
plt.gca().tick_params(axis='y', which='major', length=16, width=4)
#plt.ylim(-0.5, 1.5)
# Adding a legend
#legend = plt.legend(loc='upper right', prop={'size': 26}, bbox_to_anchor=(0.99, 0.98), borderaxespad=0.,ncol=2)
#legend.get_frame().set_linewidth(4)
#legend.get_frame().set_edgecolor('black')

ax = plt.gca()
for spine in ax.spines.values():
    spine.set_linewidth(4)

# Set x-axis limits within the data range
plt.xlim(min(dr_obs.index), max(dr_obs.index))
# Display the plot
#plt.xlim(1961, 2022)
plt.show()


# # Output

# # WGW_Max

# In[70]:


fr_train = pd.concat([fr_sim1, fr_sim2, fr_sim3, fr_sim4, fr_sim5], axis=1)
column_names = ['fr_sim1', 'fr_sim2', 'fr_sim3', 'fr_sim4', 'fr_sim5']

# Assign column names to fr_train_df
fr_train.columns = column_names

# Print the DataFrame
print(fr_train)


# In[71]:


fr_train.to_csv('fr_model.csv', index=True)


# In[72]:


fr_obs = pd.concat([fr_obs], axis=1)
obs_name = ['fr_obs']

# Assign column names to fr_train_df
fr_obs.columns = obs_name

# Print the DataFrame
print(fr_obs)


# In[73]:


fr_obs.to_csv('fr_obs.csv', index=True)


# In[74]:


#fr_train= fr_train.values.reshape(fr_train.shape[0], 1,fr_train.shape[1])


# # WGW_Min

# In[75]:


dr_train = pd.concat([dr_sim1, dr_sim2, dr_sim3, dr_sim4, dr_sim5], axis=1)


# Assuming fr_sim1 to fr_sim6 are your column names
column_names = ['dr_sim1', 'dr_sim2', 'dr_sim3', 'dr_sim4', 'dr_sim5']

# Assign column names to fr_train_df
dr_train.columns = column_names

# Print the DataFrame
print(dr_train)


# In[76]:


dr_train.to_csv('dr_model.csv', index=True)


# In[77]:


dr_obs = pd.concat([dr_obs], axis=1)
dr_name = ['dr_obs']

# Assign column names to fr_train_df
dr_obs.columns = dr_name

# Print the DataFrame
print(dr_obs)


# In[78]:


dr_obs.to_csv('dr_obs.csv', index=True)


# In[ ]:





# In[ ]:





# In[ ]:





# In[ ]:





# In[ ]:




